Electron scattering states at solid surfaces calculated with realistic potentials 
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Scattering states with low-energy-electron diffraction 
asymptotics are calculated for a general non-muffin tin po- 
tential, as, e.g., for a pseudopotential with a suitable barrier 
and image potential part. The latter applies especially to the 
case of low lying conduction bands. The wave function is de- 
scribed with a reciprocal lattice representation parallel to the 
surface and a discretization of the real space perpendicular 
to the surface. The Schrodinger equation leads to a system 
of linear one-dimensional equations. The asymptotic bound- 
ary value problem is confined via the quantum transmitting 
boundary method to a finite interval. The solutions are ob- 
tained based on a multigrid technique that yields a fast and 
reliable algorithm. The influence of the boundary conditions, 
the accuracy, and the rate of convergence with several solvers 
are discussed. The resulting charge densities are investigated. 

PACS numbers: 61.14.Hg, 2.70.Bf, 2.60.LJ, 79.60.-i 



Electron spectroscopies at low energies below 100 eV 
are sensitive to the shape of the surface near potential. 
Sophisticated bandstructure calculations involve such po- 
tentials, but almost exclusively consider bound states. 
To obtain the same accuracy for current carrying states 
is still a challenge. Here we will consider final states 
of photoemission, which are time-reversed states of low 
energy electron diffraction (LEED). Typical energies of 
interest are 10 eV and higher. With the_advent of space- 
filling full-potential LEED calculations £l the correct sur- 
face barrier can be included, but this requires, especially 
at these energies, extended CPU time. Another approach 
relies on smooth continuous matchingH of the solutions 
inside the crystal to those in the vacuum. The inclu- 
sion of ihe surface potential by the propagatijon-matrix 
methodEl has proven to be an ill-posed problem.Q The rea- 
son is the formulation as an initial value problem, which 
does not ensure the crucial continuous dependence on 
the boundary values. Tbw-,is guaranteed by the two- 
side boundary conditions B-lI Modern treatments of el- 
liptic problems often employ a discretization in direct 
space instead of choosing physically motivated basis func- 
tions. A recent work on LEED and low energy positron 
diffraction presented a three-dimensional finite-difference 



method.El For the lower energies of ultraviolet photoe- 
mission spectroscopy, the Coulomb singularities of the 
cores are less important. They can be avoided by using 
pseudopotentials, for which a mixed representation of the 
wave function seems to be suitable here. 

In this contribution we address a direct, methodically 
simple, and fast solution of the Schrodinger equation with 
the scattering asymptotics treated as introduced by Lent 
and Kirkner.Q The problem is formulated as a two-side 
boundary problem. The calculations refer exemplary to 
the GaAs (110) surface. The Schrodinger equation has 
to be solved for a given energy E and a given direction of 
the electron incident on the surface with a surface parallel 
wave vector fc|| . The translation symmetry parallel to the 
surface allows a description of the wave function $ and 
the potential V in the Laue representation, consisting of 
a Fourier decomposition in the xy plane parallel to the 
surface with the Fourier coefficients ipn and V^ depending 
on the coordinate z of the direct space perpendicular to 
the surface, 
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The two-dimensional vector p lies in the xy plane, the 
imaginary optical potential S accounts for the attenu- 
ation owing to many-particle effects and other inelastic 
losses. The summation is over the reciprocal lattice vec- 
tors g. In the Laue representation the Schrodinger equa- 
tion appears as a system of linear one-dimensional differ- 
ential equations, which we discretize by a grid with step 
size hz in the z direction. This leads to the linear ma- 
trix equation Ah^ 0h^ — with the solution iph;, and the 
coefficients 
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The indices are the reciprocal lattice vectors g, cf and 
the grid points i,i'. The complex matrix Ah^ is non- 
hermitean and indefinite. Since the potential relies on 
the Laue representation without any restrictions, the 
method allows direct application of potentials from mod- 
ern band-structure computations. Nonlocal pseudopo- 
tentials can be applied in a way similar to their use in 
three-dimensional finite difference calculations of elec- 
tronic structure.cl For the calculation here we take as 
an example a potential, which was repeatedly applied 
to photoemission calculations of GaAs(110).Ll It is a lo- 
cal pseudopotential with the surface barrier taking into 
account relaxation, corrugation, and a smooth saturated 
image potential. 

Equation ^ requires boundary conditions at the bulk 
and at the vacuum side of the grid. The simplest model 
of a surface potential is a single step towards the vac- 
uum, as used in matching calculations. A first guess for 
boundary conditions far away from the surface inside the 
crystal and that outside within the vacuum are Dirich- 
let values taken from such preliminary calculations. In 
Fig. n the modulus of the wave functions for different 
grid sizes is plotted. Plots (c) and (d) show irregularities 
around the bulk boundary. Even with Neumann or mixed 
boundaries this feature cannot be suppressed. Because 
the damping reduces the wave functions to zero deep in 
the bulk, the required asymptotic behavior is fulfilled by 
taking a zero value at a boundary sufficiently remote from 
the surface. Since the damping depends on energy, the 
coordinate of the boundary in the bulk should be cho- 
sen as energy dependent and automatically adjusted by 
tracking the neighboring values of the wave function. 

In the vacuum the wave functions still differ even for 
high grid sizes as shown in Fig. |l| (a) and (b). Therefore, 
we investigate this boundary more closely. In a LEED 
experiment a beam of electrons is scattered at the sur- 
face of the crystal. The vacuum boundary condition has 
to guarantee one single incoming plane wave. This can be 
tested by a half-sided Fourier transformation of the wave 
function in the vacuum, assumin that the Fourier grid is 
well separated from the vicinity of the surface. The solu- 
tion from a matching calculation does fulfill the incom- 
ing beam asymptotics, the subsequent grid calculation 
with the same step potential does not. The Fourier spec- 
trum of the grid wave function shows additional incom- 
ing ghost waves that have contributions of up to 15% and 
thereby strongly deteriorate the physical situation. Neu- 
mann and mixed boundaries from the matching do not 
improve this unsatisfactory situation. Again the match- 
ing does not provide a correct boundary condition even 
far away from the crucial vicinity of the surface. 

The correct boundary conditions in the vacuum are 
achieved by use of the elegant and simple quantum trans- 
mitting boundary method (QTBM)J3 The wave function 
in the vacuum is expanded into a set of propagating waves 
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.^ . and (jy^ are the amplitudes of the given incom- 
ing and unknown outgoing waves, respectively. The 
grid wave function '^ grid is that of Eq. (P. For a 
short description of the quantum transmitting bound- 
ary method let LN{ ) be a linear function and FT{ ) the 
two-dimensional Fourier transformation. Continuity at 
boundary z,. then yields at that plane 
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From the continuity of the normal derivative follows 
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This mixed linear boundary value problem is the QTBM 
that is inserted in Eq. (^ with the derivative being dis- 
cretized. The implementation of the QTBM is simple 
here due to the Laue representation. With the QTBM 
the calculated wave function fulfills the correct LEED 
asymptotics in the case of the step barrier as it does for 
any corrugated barrier too. Figure y shows the Fourier 
spectrum normal to the surface of a grid calculation with 
the QTBM basing either on a step potential or basing on 
the true barrier. The negative wave numbers correspond 
to incoming waves. Both solutions show the postulated 
LEED asymptotics of one incoming wave. The barrier 
potential influences strongly the solution, which is also 
illustrated by the plot of the wave functions in Fig. 0. 

Equation (|^) with the boundary conditions inserted 
gives a quadratic, linear, and inhomogeneous system of 
equations, which has a block tridiagonal coefficient ma- 
trix. As an example we consider the GaAs(llO) surface 
with a normal incident electron beam at an energy of 
Ef = 18 eV, which in vacuum corresponds to a kinetic 
energy of Ekin = 12.75 eV. The lateral plane wave ba- 
sis of 57 reciprocal lattice vectors proved to be sufficient 
for the employed potential. With the lattice constant of 
a = 5.654 A, the boundaries are chosen to be at ±15a, 
enclosing 45 layers with 90 atoms. The large distance 
from the surface to the vacuum boundary ensures a suf- 
ficient decay of the image potential. The positions of the 
boundaries deep inside the crystal and far outside in the 
vacuum cause a large grid, whose step size h^ is bounded 
by several conditions. First, the wave function has to 
be correctly represented on the grid. This bound can be 
estimated by use of the sampling theorem. Furthermore, 
a small step size is needed for a stable discretization and 
a safe convergence of the iterative solvers. The Laplace 
operator has to dominate the zero-order terms, which 
require h^ < 0.06a here. Further bounds are given by 
the demands of applications. With h^ = 0.012a a typical 
number of 2500 grid points and 142 500 equations results. 

Since the coefficient matrix results from an elliptic 
equation, is sparse, and is additionally a band matrix, 
a variety of solvers are available. As for the photoemis- 



sion spectra a large series of final states has to be calcu- 
lated, it is very important to reduce the CPU time per 
final state. Several direct and iterative solvers have been 
tested. Some of the resulting CPU time and memory 
requirements are given in Table I. Time and memory in- 
crease linear with the number of grid points due to the 
simple structure of the coefficient matrix. The fastest 
direct solver was a routine in band storage mode per- 
forming an LU factorization that needs 98 s for the test 
problem with a slope of 3.9 x 10~^s per grid point. It 
needs less than half of the time than a LEED calculation 
with space-filling potentials at this energy or a conven- 
tional matching at a potential step. 

Iterative solvers in combination with multigrid 
methods lase highly successful in solving differential 



equations.LL!J They have recently been applied, to the cal- 
culations of large-scale electronic structure.Eil We imple- 
mented a two- and a three-grid method. For the smooth- 
ing, the squared Jacobi iteration has been applied. With 
a damping factor of 0.5 the best convergence rates have 
been obtained. The Laue representation simplifies the 
implementation of the QTBM, but it causes strong co- 
diagonals due to the potential coefficients. Therefore on 
the first level the equations must be treated by a direct 
solver, and the Jacobi iteration is only used for smoothing 
on the finer grids. For the step size of 0.018a on the coarse 
grid the two-grid iteration needs 20 steps to achieve the 
same mean defect of '-^ 10^^^ as the direct solver. The 
three-grid methods do not improve the convergence rates. 
The defect after 20 iterations was for a W cycle 6xl0~^° 
and for a V cycle 5 x 10~^. Theoretically the multigrid 
method needs 0{N) operations to solve a system of N 
linear equations. The direct method for band matrices 
needs 0{Nw^) operations. Because the number of co- 
diagonals w is the constant number of reciprocal lateral 
lattice vectors, both methods have the same asymptotic 
behavior. However, if a lower accuracy of the solution 
is acceptable, the multigrid method becomes more favor- 
able. If, e.g., a defect of 10~^ is sufficient, the two-grid 
method needs five iterations and 91 s and gives an error 
in the modulus of the wave function of less than 0.2%. 
Additionally, the maximum memory used can be reduced 
within the multigrid method. If the system of equations 
is solved directly, a (Sw + l)xN matrix containing the 
LU factorization on the fine grid is needed. In the multi- 
grid iteration with depth I the (2-u;-t-l)xiV matrix on the 
fine grid and the (Siu+l) x (7V2~('~^^) LU matrix on the 
coarse grid are stored. From the three-grid iteration on- 
ward, the required memory for the multigrid procedures 
becomes less than in the other methods. The multigrid 
calculations may possibly be accelerated by the use of 
iterative solvers better suited to the indefinite problem 
than the squared Jacobi iteration. 

As a first application the charge-density distribution 
of a LEED state is shown in Fig. 0. The interface crys- 
tal vacuum and the exponential decay of the wave func- 
tion into the crystal are clearly visible. In the vacuum 
a distance of several lattice constants from the surface 



is necessary to evolve an interference pattern typical of 
superimposed waves, c.f. Fig. g. The importance of the 
correct treatment of the potential barrier, which is al- 
ready obvious from the different LEED intensities in Fig. 
0, is illustrated by the strong charge fluctuations in the 
surface region. In contrast to tracking multiple scattering 
paths in conventional LEED calculations, the properties 
of LEED states are much easier interpreted with a single 
wave function. It is interesting to see that there are well 
localized regions of high charge density even at this non- 
bonding energy. In applications to photoemission these 
regions will give strong contributions to the photocur- 
rent. Thus, this treatment supports a local interpreta- 
tion of the emission. First photoemission calculations 
with these final states were promising. 

Solutions of the Schrodinger equation with scattering 
boundary conditions have been calculated with a realistic 
potential of the GaAs (110) surface. The correct asymp- 
totic behavior was obtained by the quantum transmitting 
boundary method. The Laue representation allows a sim- 
ple implementation of the boundary conditions. Several 
multigrid methods and direct solvers have been tested. 
The fastest were a routine in band storage mode with 
LU factorization and, when stopping at a slightly lower 
accuracy, a two-grid method. The multigrid method is 
competitive with direct solvers due to its higher flexibil- 
ity. A general implementation of the potential allows us 
to use potentials from modern ab-initio techniques. In 
an application, the calculated charge density of a LEED 
state showed distributions with localized accumulation 
regions of the excited electron. It may yield a comfortable 
access to a local interpretation of photoemission spectra. 
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Method 


CPU time (s) 


Memory (MB) 


Sparse Gauss 


1862 


747 


Sparse LU 


(1862) 


>2000 


Band LU with refinement 


257 


654 


Band LU without refinement 


98 


392 


Two grid 


208 


458 


Three grid V cycle 


343 


360 


Three grid W cycle 


270 


360 




TABLE I. CPU time and memory used by different di- 
rect and iterative routines for a system of 142 500 equations 
on a CRAY-J916. As direct solvers routines from the Inter- 
national Mathematics and Scientific Library were used. For 
the routine with the LU factorization of a sparse matrix the 
time is extrapolated from tests with smaller grids. The given 
memory is required by the 64-bit architecture of the Cray. 



FIG. 1. Sensitivity of the solution to the range of cal- 
culation: wave functions at {x, y) — (0, 0) for different grid 
sizes. The boundary positions are marked in a, at which 
boundary conditions from a previous matching calculation are 
used. The potential for the matching was a step, here it is 
the smooth barrier. 



1.0 



0.2 



0.1 



0.0 



^ 



I potential barrier 
Q potential step 



y1 



V 



„jUL. 



-3.0 



-2.0 -1.0 0.0 1.0 2.0 

z - cornponent of the wave vector (2ii/a) 



3.0 



FIG. 2. Comparison of grid calculations with a step po- 
tential and a potential barrier: modulus of the Fourier ampli- 
tudes normal to the surface in vacuum at (a;, y) = (0, 0). The 
squared amplitudes give the LEED-intensities. 
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FIG. 3. Comparison of the wave functions from a step po- 
tential and a potential barrier at {x,y) — (0,0). The position 
of the step is indicated by the arrow. 




FIG. 4. Gray scale image and contour plot of the charge 
density in the yz plane at a:: = with a z interval from — 3.55o 
to 3.95a. In the y direction two unit cells are shown. The 
atoms lying in the plane are drawn black, and in the right 
cell the projected positions of the remaining atoms are gray. 
The different scales for the y and z directions are indicated 
by the Angstr0m scale. 



